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Time-domain  solution  of  higher-order  parabolic 
equations 

Michael  D.  Collins 

Naval  Ocean  Research  and  Development  Activity 
Stennis  Space  Center,  MS  39529 


Abstract 


A  higher-order  time-domain  parabolic  equation  (TDPE)  is  derived  from  a  Fade  series[7],  solved 
numerically,  and  applied  to  underwater  acoustics  problems.  The  higher-order  TDPE  solution  is 
accurate  for  problems  involving  very  wide  propagation  angles  and  large  variations  in  sound  speed. 
Its  applications  include  propagation  near  the  source,  propagation  out  to  very  long  ranges,  and 
propagation  over  a  hard  ocean  bottom.  The  higher-order  TDPE  is  valid  in  both  shallow  and  deep 
water.  The  accuracy  of  the  model  is  demonstrated  with  benchmark  calculations.  The  model  is 
applied  to  illustrate  mode  cutoff  in  a  range-dependent  ocean. 

1.  Introduction 


The  parabolic  equation  (PE)[1]  is  a  verj  -cful  model  for  range-dependent  propagation  calcu¬ 
lations  in  underwater  acoustics.  The  PE  was  originally  accurate  only  for  problems  involving  limited 
bottom  interaction  and  narrow  propagation  angles.  However,  the  PE  was  later  extended  to  handle 
wide-angle  propagation  [2-4]  and  bottom  interaction  [3-5].  Benchmark  comparisons  showed  that 
the  wide-angle  PE  is  accurate  for  many  underwater  acoustics  problems  [6],  The  accuracy  of  the 
PE  has  recently  been  further  improved.  A  higher-order  PE  based  on  a  Pade  series  [7,8]  produces 
solutions  as  accurate  as  the  outgoing  coupled  mode  solution[9]  for  very  wide  propagation  angles 
and  large  differences  in  sound  speed.  In  particular,  the  higher-order  PE  is  accurate  for  propagation 
very  near  the  source,  propagation  out  to  very  long  ranges,  and  propagation  over  high-speed  ocean 
bottoms. 

Two  approaches  exist  for  applying  the  parabolic  approximation  to  pulse  propagation.  The  PE 
has  been  applied  to  solve  time-domain  problems  in  the  frequency  domain  using  Fourier  synthesis. 
[10,11]  This  is  probably  not  an  optimal  approach  if  the  solution  is  desired  at  many  points  in  the 
domain.  If  a  sequence  of  snapshots  of  the  solution  were  desired,  for  example,  one  would  have  to 
perform  Fourier  synthesis  at  each  point  in  each  snapshot.  The  time-domain  parabolic  equation 
(TDPE),  which  is  the  inverse  Fourier  transform  of  the  PE,  has  been  developed  for  range-dependent 
time-domain  propagation  calculations.[l2-16]  With  this  approach,  one  automatically  obtains  the 
solution  at  all  points  in  the  domain. 

In  this  paper,  a  higher-order  TDPE  is  derived,  solved  numerically,  and  applied  to  underwater 
acoustic  pul?«'  ropagation.  In  addition  to  allowing  one  to  avoid  Fourier  synthesis,  the  higher- 
order  TDPE  has  all  of  the  advantages  of  the  higher-order  PE.  The  numerical  solution  is  based 
on  Galerkin’s  method  because  the  depth  operators  are  relatively  complicated.  The  higher-order 
TDPE,  which  is  accurate  for  propagation  in  both  shallow  and  deep  water,  is  compared  with  a  wide- 
angle  TDPE  designed  for  shallow  water.  The  accuracy  of  the  higher-order  TDPE  is  demonstrated 
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with  benchmark  comparisons.  The  model  is  used  to  illustrate  mode  cutoff  in  an  ocean  with  an 
upward  sloping  bottom. 

2.  Factoring  the  wave  equation 

We  work  in  cylindrical  coordinates  with  r  being  the  range  from  a  point  source  and  z  being  the 
depth  below  the  ocean  surface.  We  assume  for  now  that  the  domain  is  stratified  and  that  azimuth 
dependence  can  be  neglected.  We  remove  the  spreading  factor  r~ 5  from  the  acoustic  pressure  p 
and  begin  with  the  farfield  wave  equation 

&P  &P  _  J_f^P  (2  11 

dr 2  dz2  p  dz  dz  c 2  dt2  ' 

where  t  is  time,  c  is  sound  speed,  and  p  is  density.  Equation  (2.1)  can  be  factored  in  two  ways: 


0  The  TDPE,  which  is  the  inverse  Fourier  transform  of  the  PE,  is  obtained  by  approximating  the 

square  root  in  Eq.  (2.2).  The  progressive  wave  equation  (PWE)  is  derived  from  Eq.  (2.3).  The 
TDPE  and  PWE  are  valid  to  leading-order  in  domains  in  which  range  dependence  is  a  perturbation. 

A  numerical  solution  of  the  PWE  based  on  the  method  of  alternating  directions  and  nonlinear 
capability  for  the  PWE  have  been  developed,  [14]  and  the  PWE  has  been  extended  to  handle 
density  variations  and  sediment  attenuation. [15]  However,  the  PWE  can  not  be  extended  to  handle 
wide-angle  propagation,  and  the  TDPE  is  easier  to  initialize  and  is  better  suited  to  handle  range- 
dependent  problems. [16]  Thus  recent  development  has  involved  the  TDPE,  which  has  been  extended 
to  handle  sediment  dispersion  and  wide-angle  propagation  in  shallow  water.[16] 

3.  The  higher-order  TDPE 


The  1-term  Taylor  series 
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Table  I:  Comparison  of  Taylor  and  Fade  series. 


X 

4- term 
Taylor 

1-tcrm 

Pade 

2- term 

Pade 

3-term 

Pade 

■/l  +  x 

0.25 

1.11801 

1.11765 

1.11803 

1.11803 

1.11803 

0.50 

1.22412 

1.22222 

1.22472 

1.22474 

1.22474 

0.75 

1.31870 

1.31579 

1.32274 

1.32287 

1.3228a 

1.00 

1.39844 

1.40000 

1.41379 

1.41420 

1.41421 

1.25 

1.45639 

1.47619 

1.49904 

1.49996 

1.50000 

1.50 

1.48193 

1.54545 

1.57931 

1.58105 

1.58114 

1.75 

1.46078 

1.60870 

1.65523 

1.65812 

1.65831 

2.00 

1.37500 

1.72000 

1.79584 

1.80221 

1.80278 

2.50 

0.91943 

1.76923 

1.86124 

1.86994 

1.87083 

2.75 

0.49545 

1.81481 

1.92376 

1.93519 

1.93649 

3.00 

-0.10156 

1.85714 

1.98361 

1.99817 

2.00000 

has  been  used  to  derive  the  narrow-angle  TDPE  and  the  narrow-angle  PWE,  which  are  accurate 
for  propagation  angles  up  to  about  15  degrees.  The  1-term  Pade  series 

_ _ _  I, 

\ZT+i  -  1  =  — 2-r-  +  0(x3)  (3.2) 

1  +  \x 

has  been  used  to  derive  a  wide-angle  TDPE  for  shallow  water,  which  we  refer  to  as  TDPES  and 
which  is  accurate  for  propagation  angles  up  to  about  40  degrees.  To  obtain  a  higher-order  TDPE 
that  is  accurate  for  the  very  wide  propagation  angles  that  occur  near  the  source  and  over  a  hard 
ocean  bottom,  an  approximation  for  the  square  root  function  that  is  very  accurate  for  x  2  1  is 
required.  Since  the  Taylor  series  converges  only  for  |  x  |<  1,  many  terms  are  required  for  x  =S  1. 

The  following  generalization  of  Eq.  (3.2)  is  not  restricted  to  |  x  |<  1 


-1  =  E  rzr^  +  o(x^), 

where  n  is  the  number  of  terms  in  the  Pade  series  and 

2  .  2  i* 

n  -  —  _  cm*  _ _ 


2  71  -f-  1  2 71  -f  1 

f  \ 

2 

(» 

—  cos2  - . 

2n  +  1 

(3.3) 


(3.4) 

(3.5) 


Availability  Codes 


iDlat 


Avail  and/or 
Special 


□  □ 


Normalized  pressure  Normalized  pressure  Normalized  pressure 
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Relative  timefms) 


Figure  1:  Time  series  at  r  =  500m  and  z  =  200m  for  a  Gaussian  pulse  in  a  waveguide  with 
perfectly  reflecting  boundaries.  The  dashed  curves  are  the  image  solution.  The  solid  curves  are  the 
(aj  TDPEi,  (b)  TDPE2,  and  (c)  TDPE3  solutions. 
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Since  the  Pade  series  is  valid  outside  the  radius  of  convergence  of  the  Taylor  series,  relatively  lew 
terms  are  needed  for  x  =  1.  We  illustrate  this  in  Table  I.  The  4-term  TayJor  series  is  better  than 
the  1-term  Pade  series  for  i  <  1,  but  the  1-term  Pade  series  is  better  for  x  >  1.  The  2-term  Pade 
series  and  the  4-term  Taylor  series  are  both  correct  to  0( xs)  for  small  x.  Yet  the  2-terin  Pade  series 
is  substantially  better  than  the  4-term  Taylor  series.  The  3-term  Pade  series  is  fairly  accurate  well 
beyond  the  radius  of  convergence  of  the  Taylor  series  near  i  =  3. 

The  effects  of  attenuation  and  dispersion  are  less  important  than  refraction  and  diffraction 
for  most  problems.  Thus  we  do  not  derive  corrections  for  the  attenuation/dispersion  operator  in 
TDPE,.  For  now,  we  assume  that  c  is  real  (no  attenuation)  and  independent  of  u  (no  dispersion) 
in  the  analysis.  From  Eqs.  (2.2)  and  (3.3),  we  obtain  the  following  higher-order  TDPE 


where 


Oj.n  —  ttj .  H  • 

We  define  u(r,z,t)  =  p(r,z,t  +  r/c o)  and  rearrange  Eq.  (3.6)  to  obtain 


(3.6) 


(3.t: 


(3.ft) 


(3.!)) 


(3.10) 


(3.1!) 


We  refer  to  Eq.  (3.11)  as  TPDE„.  In  the  derivation  of  TDPEj,  it  was  assumed  that 
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/  1  1  \  d2u  d2u  1  dpdu 

\c2  Cq)  dt2  <<:  dz2  pdzdz 


(3.12) 


which  is  valid  in  shallow  water.  Since  this  assumption  was  not  made  to  derive  TDPE„,  we  deduce 
that  TDPEi  is  the  generalization  of  the  wide-angle  PE  to  deep  water. 


4.  Numerical  solution 


The  alternating  directions  solution  of  TDPEn  requires  numerical  solutions  for  each  of  the 
following  n  +  1  equations 


JFl  ,  s  (—  -  = 

3'n  drdt  J,n  \dz2  pdzdzj  dr 


d2 

a? 


i  dp 

POz 


du 


('t.i) 


(4.2) 


Since  Eq.  (4.1)  is  similar  to  the  refraction  term  of  TPDES,  and  Eq.  (4.2)  is  similar  to  the  diffraction 
term  of  TDPES,  the  numerical  solution  developed  in  Ref.  16  can  be  modified  slightly  to  obtain  the 
numerical  solution  of  Eq.  (3.11).  Without  the  rearrangement  of  Eq.  (3.6),  it  would  be  necessary 
to  solve  n  equations  similar  to  Eq.  (4.2)  as  well  as  n  third-order  equations  that  are  much  more 
complicated  than  Eq.  (4.1). 

As  in  Ref.  16,  the  source  function  f(t)  is  assumed  to  have  compact  support,  and  a  time  window 
ti  <  t  <  t2  that  contains  the  signal  at  all  times  is  chosen.  The  boundary  condition  u  =  0  is  imposed 
at  the  pressure  release  surface,  deep  within  the  sediment  at  z  —  z,\i  from  which  no  energy  returns 
to  the  water  column  due  to  attenuation,  and  after  the  signal  has  passed  the  obsever  at  t  =  <2-  The 
boundary  conditions  u  =  du/dl  =  0  are  imposed  before  the  signal  is  detected  at  t  =  t\.  Equation 
(4.1)  is  a  first-order  hyperbolic  equation  that  can  be  solved  with  the  Lax-Wendroff  scheme[l7]. 

Galerkin’s  method  is  used  to  discretize  depth  dependence  in  Eq.  (4.2).  The  resulting  equation 
is  then  solved  with  Crank-Nicolson  integration  in  r  using  centered  differences  in  t  while  sweeping 

from  t  =  t [  to  l  =  <2-  We  define  the  depth  grid  points  a,  =  lAz.  The  basis  functions  4*,(^)  vanish 
for  |  z  —  Zj  |>  A z,  increase  linearly  from  0  to  1  over  z,_i  <  z  <  z,,  and  decrease  from  1  to  0  over 
Zi  <  z  <  Zj+i.  We  define  u,(r,t)  =  u(r,zitl)  as  well  as  Qj  =  0(zj)  and  3>,  =  ${z,)  for  arbitrary 
functions  0  and  <P.  The  basis  functions  provide  the  approximations 
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u(r,z,t)  2  E  «i(r>0*i(*)  (4-1) 

t 

*(*)  a  E  (-«--*) 

i 

e(z)  a  E  ««*.-(*)•  ('‘-5) 


The  depth  operator  Qz  is  discretized  with  Galerkin’s  method,  as  follows: 


Qz*  U 


f*,Qz*dz 

JViJz 


(-1.6) 


Substituting  Eqs.  (4.3),  (4.4),  and  (4.5)  into  Eq.  (4.6),  we  obtain  the  following  approximations  for 
the  depth  operators: 


Qu  u=2,  a  Q-11./  + 

0,_1  +  60j  +  0;+1  ©,+  l  +  0, 

- - -  tt,  +  - L-LJ2 -  + 


(4.7) 


0 


d2u 

8? 


0, 

(A;)2 


20,  .  0, 

5  +  T-  To  Mi+1 


(AO2  ( Az)2 


(4.8) 


*£&$!(  I  a,  (&-i  +  2»0  (0i-i  -  0,)  + 

'*=*■  -  6(^7  “-1  + 

^i-i(9j  -  Qj-i)  4  2$, (20,  -  0,-i  -  0,+i )  4-  ^i+i(0j  -  0,4-1 ) 

6(Al? 

(^■+i  +  24>,)  (0,+i  -  0,) 

6(a7?  u,+'- 


(4.9) 


5.  Examples 

To  demonstrate  the  ability  of  TDPE„  to  handle  very-wide-angle  propagation,  we  consider 
a  waveguide  of  thickness  300m  with  pressure-release  top  and  bottom  boundaries  in  which  <•  = 
1500m/s.  The  Gaussian  source  /(<)  =  exp  [-(i/t)2]  is  placed  at  z  =  25m,  where  //  =  150s_l.  The 
image  solution,  which  is  exact,  is  used  to  initialize  the  field  at  r  =  200m,  and  we  take  cq  —  1500m/s. 
The  TDPEi,  TDPE2,  and  TDPE3  solutions  are  compared  with  the  image  solution  in  Figure  1. 

Each  of  the  solutions  is  very  accurate  for  the  first  arrivals,  which  propagate  at  small  angles. 
However,  the  agreement  improves  with  n  for  the  later  arrivals,  which  propagate  at  larger  angles. 
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Range  (m)  —  2000.00  Range  (m)  =  4000.00 


Range  (m)  =  6000.00 


Range  (m)  =  8000.00 


Range  (m)  =  10000.00 
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S00  400  300  800  100 

Relative  time  (ms) 


Figure  2:  Contour  plots  of  a  Gaussian  pulse  in  a  refracting  ocean. 
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2hl> 


Range  (km) 


Figure  3:  Transmission  loss  at  z  =  390m  for  a  50Hz  source  in  a  refracting  ocean.  The  dashed  curves 
are  the  wide-angle  PE  solution.  The  solid  curves  are  the  (a)  narrow-angle  PE,  (b)  TDPEj,  and  (c) 
TDPEi  solutions. 
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Range  (m)  =  5000.00 


Range  (m)  =  7000.00 


Figure  A:  Contour  plots  of  a  Manning  weighted  sinusoidal  pulse  in  a  range-dependent  ocean 
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In  past  studios  of  the  TDPE,  a  stability  condition  for  the  numerical  solution  of  the  refraction 
operator  has  been  discussed.  However,  the  numerical  solution  of  the  diffraction  operator  appeared 
to  be  unconditionally  stable  based  on  numerical  results.  While  performing  the  calculations  for  the 
previous  example,  however,  we  discover.:.!  a  new  stability  condition  involving  the  grid  spacings  A z 
and  At.  In  a  homogeneous  medium,  the  numerical  solution  of  the  diffraction  operator  is  unstable 
for  A^  =  coAt  and  n  >  1.  The  solution  appears  to  be  stable  for  all  n  if  A z  t>  AcoAt,  where 
numerical  experiments  give  A  =  1.4. 

To  demonstrate  the  ability  of  TDPE„  to  handle  large  variations  in  sound  speed,  we  consider 
an  ocean  of  depth  400m  in  which  c  increases  linearly  from  1500m/s  at  z  —  0  to  1600  m/s  at  z 
=  100  m  (this  is  not  intended  to  represent  a  realistic  sound-speed  profile).  In  the  sediment,  c  = 
I700m/s,  p  =  1  .5g/cm3,  and  the  attenuation  is  0  —  0.5dR/A.  The  Gaussian  source  function  with 
v  —  150s-1  is  placed  at  z  —  50m,  and  we  take  Co  =  1500m/s.  The  homogeneous  half-space  field 
[18]  is  used  as  an  initial  condition  at  r  =  100m.  The  plane-wave  loss  operator  of  Ref.  16  is  used 
to  model  attenuation.  However,  we  have  found  that  greater  accuracy  is  obtained  by  using  c  rather 
than  c0  in  the  loss  operator.  Sediment  dispersion  is  neglected. 

A  sequence  of  contoin  plots  of  p  computed  with  TDPEi  appears  in  Figure  2.  Solid  contours 
correspond  to  p  >  0;  dashed  contours  correspond  to  p  <  0.  The  solid  line  marks  the  ocean  bottom. 
The  response  to  /  is  convolved  as  in  Ref.  16  to  obtain  the  response  to  a  50Hz  time-harmonic  source. 
Transmission  loss  for  the  iDPEj,  Tl)PEa.  wide-angle  PE,  and  narrow-angle  PE  solutions  appears 
in  Figure  3.  The  narrow-angle  PE  solution  has  large  phase  errors.  The  TDPE3  solution  is  better, 
but  it  too  has  a  large  error  due  to  the  strong  refraction.  The  excellent  agreement  of  the  TDPEi  and 
wide-angle  PE  solutions  demonstrates  the  ability  of  TDPEi  to  accurately  handle  pulse  propagation 
in  deep  water  and  shows  that  the  plane-wave  loss  operator  is  accurate  for  this  problem. 

The  phenomenon  of  mode  cutoff  in  a  range-dependent  oman  has  been  illustrated  for  time- 
harmonic  signals. [19]  To  illustrate  energy  cutoff  in  the  time  domain,  we  apply  TDPE]  in  a  range- 
dependent  ocean  in  which  the  depth  (I  is  defined  by 


il(  r  I 


100m  for  r  <  5km 

(l  -  r )  100m  for  5km  <  r  <  10km 
50m  for  r  >  1 0km. 


(6.1) 


We  take  c  =  1500m/s  in  the  water  and  c  =  1600m/s,  p  —  i.5g/cm3.  and  fi  —  0.5dB/A  in  the 
sediment.  The  source  function  is  the  Hanning  weighted  sinusoid  [20] 


((l  -  cos  f  uit\  sin  off  0  <  t  <  T 

(6.2) 

0  otherwise 

where u;  =  1007T.S-1  and  T  =  H-t/uj.  Since  this  source  function  has  a  fairly  narrow  frequency  content, 
it  is  effective  for  illustrating  normal  mode  behavior.  The  source  is  placed  at  -  =  25m.  and  we  apply 
the  half-space  field  at  r  =  100m. 
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Snapshots  of  the  TDPEi  solution  appear  in  Figure  4.  The  first  three  modes  arc  dearly  visible 
at  r  =  5kin.  The  vertical  dependence  of  the  field  matches  the  first  mode  for  the  first  arrivals  as 
only  one  extremum  appears.  A  maximum  and  a  minimum  occur  near  t  =  150ms  due  to  arrivals 
of  the  second  mode.  The  third  mode,  which  has  three  extrema,  is  apparent  for  t  >  200ms.  Cutoff 
occurs  for  two  of  the  modes  as  d  decreases.  The  third  mode  is  no  longer  evident  at  r  =  9km.  The 
second  mode  disappears  by  r  =  1 1km. 

6.  Conclusion 


A  higher-order  TDPE  has  been  derived  and  solved  numerically.  Since  the  model  is  based  on 
a  Pade  series  that  is  very  accurate  and  is  not  restricted  to  the  unit  circle  in  the  complex  plane, 
it  is  accurate  for  problems  involving  very  wide-angle  propagation  as  well  as  large  sound-speed 
variations.  Since  the  higher-order  TDPE  splits  into  terms  similar  to  the  refraction  and  diffraction 
terms  of  TDPE,,  it  can  be  solved  with  similar  methods.  Although  the  coefficients  of  the  terms  in 
the  higher-order  TDPE  are  relatively  complicated,  they  are  easily  handled  with  Galerkin’s  method. 
The  higher-order  TDPE  solution  is  comparable  in  accuracy  to  outgoing  coupled  normal  mode 
solutions.  In  particular,  the  model  is  valid  in  deep  water. 
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